Drivers of variation in the community response to temperature change across realms

(using mixed effects models)

Collaborators: Shane Blowes, Jon Chase, Helmut Hillebrand, Michael Burrows, Amanda Bates, Uli Brose, Benoit Gauzens, Laura Antao Assistance: Katherine Lew, Josef Hauser

Introduction

• Climate change is driving a widespread reorganization of ecological communities around the world (Parmsesan & Yohe 2003, Poloczanska et al. 2013),
• but the impacts of climate change vary substantially from one location to another and among taxa (Molinos et al. 2016 NCC, Antao et al. 2020 NEE).
• Community reorganization is substantially more common than an aggregate loss or gain of species (Dornelas et al. 2014 Science, Blowes et al. 2019 Science, Hillebrand et al. 2017 J Appl Ecol)
• There are many hypotheses for why some communities are more sensitive to warming than others, including differences in
  • metabolic rates (Dillon et al. 2010 Nature),
  • thermal physiology (Deutsch et al. 2008 PNAS, Pinsky et al. 2019 Nature),
  • microclimate availability (Burrows et al. 2019 NCC, Suggitt et al. 2018 NCC),
  • species mobility (Poloczanska et al. 2013 NCC, Burrows et al. 2011 Science, Sunday et al. 2012 NCC)
  • or generation time (Beaugrand et al. 2009 DSR II, Poloczanska et al. 2013 NCC),
  • consumers vs. producers (Petchey et al. 1999 Nature)
  • community composition (Stuart-Smith et al. 2015 Nature, Beaugrand et a. 2015 NCC, Trisos et al. 2020 Nature),
  • ecosystem productivity (Thomas et al. 2017 GCB, Brett 1971 Am Zoo),
  • exposure to human impacts (White & Kerr 2006 Ecography)
  • and among realms (Antao et al. 2020 NEE).
• Scaling up from organismal effects to whole ecological communities is complex, and yet these scales are critical for ecosystem functioning and human well-being.
• There is a need for a comprehensive test to understand where warming is driving and is likely to drive the most dramatic community turnover

Methods

• BioTime dataset, gridded to 96 km2 hexagons, summarized as temporal turnover (Blowes)
  • Temporal slope of Jaccard turnover compared to the first year
  • Same for Jaccard total
  • and Morisita-Horn turnover
• Tested explanatory variables for differences in rate of turnover:
  • Temperature trend over the time-frame of each time-series (CRU TS 4.03 on land and in freshwater, ERSST v5 in the ocean)
  • Seasonality as a metric of thermal sensitivity (Deutsch et al. 2008 PNAS). Standard deviation of monthly temperatures.
  • Average temperature as a metric of metabolic rates (Dillon et al. 2010 Nature, Antao et al. 2020 Nat E&E)
  • Mobility calculated from body mass and taxonomic group classifications of mobility mode (fly, run, swim, crawl, sessile). Fly/run/swim followed the allometric relationship in Hirt et al. 2017 Nat E&E. Crawl set at 0.1 km/hr, sessile set to 0 km/hr. Then calculated averaged within each assemblage.
  • Net primary productivity (NPP) from the merged land/ocean product produced by the Ocean Productivity group at Oregon State using methods from Zhao et al. 2005 and Behrenfeld & Falkowski 1997.
• TO DO:
  • Generation time calculated from body mass and endotherm vs. ectotherm classifications, following McCoy & Gillooly 2008 ELE. Averaged across species within each assemblage.
  • Human impact calculated from Carsten Meyer’s ecosystem cube data
  • Thermal bias calculated from Species Temperature Indices (Mike Burrows)
  • Microclimates calculated from WorldClim and BioOracle (Laura Antao)
  • Consumer vs. producer classification
• Maybe to do:
  • Species pool richness
• Differences in temporal turnover (response variable) modeled with a linear mixed effects model (nlme package, lme() function). See below for details.

library(data.table) # for handling large datasets
library(ggplot2) # for some plotting
#library(lme4)
library(nlme) # for ME models
library(beanplot) # for beanplots
library(maps) # for map
library(ggeffects) # marginal effect plots
library(gridExtra) # to combine ggplots together
library(grid) # to combine ggplots together
library(gridExtra)

options(width=500) # turn off most text wrapping

# tell RStudio to use project root directory as the root for this notebook. Needed since we are storing code in a separate directory.
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file()) 

# Turnover and covariates assembled by turnover_vs_temperature_prep.Rmd
trends <- fread('output/turnover_w_covariates.csv.gz')

# set realm order
trends[, REALM := factor(REALM, levels = c('Freshwater', 'Marine', 'Terrestrial'), ordered = FALSE)]

# group Marine invertebrates/plants in with All
trends[, taxa_mod2 := taxa_mod]
trends[taxa_mod == 'Marine invertebrates/plants', taxa_mod2 := 'All']

Log-transform some variables, then center and scale.

trends[, tempave.sc := scale(tempave)]
trends[, tempave_metab.sc := scale(tempave_metab)]
trends[, seas.sc := scale(seas)]
trends[, microclim.sc := scale(log(microclim))]
trends[, temptrend.sc := scale(temptrend)]
trends[, temptrend_abs.sc := scale(log(abs(temptrend)))]
trends[, npp.sc := scale(log(npp))]
trends[, mass.sc := scale(log(mass_mean_weight))]
trends[, speed.sc := scale(log(speed_mean_weight+1))]
trends[, lifespan.sc := scale(log(lifespan_mean_weight))]
trends[, thermal_bias.sc := scale(thermal_bias)]
trends[, consumerfrac.sc := scale(consfrac)]
trends[, endothermfrac.sc := scale(endofrac)]
trends[, nspp.sc := scale(log(Nspp))]
trends[, human.sc := scale(human)]

Do the variables look ok?

# histograms to examine
cexmain = 0.6
par(mfrow = c(3,5))
invisible(trends[, hist(tempave.sc, main = 'Environmental temperature (°C)', cex.main = cexmain)])
invisible(trends[, hist(tempave_metab.sc, main = 'Metabolic temperature (°C)', cex.main = cexmain)])

invisible(trends[, hist(seas.sc, main = 'Seasonality (°C)', cex.main = cexmain)])
invisible(trends[, hist(microclim.sc, main = 'log Microclimates (°C)', cex.main = cexmain)])

invisible(trends[, hist(temptrend.sc, main = 'Temperature trend (°C/yr)', cex.main = cexmain)])
invisible(trends[, hist(temptrend_abs.sc, main = 'log abs(Temperature trend) (°C/yr)', cex.main = cexmain)])

invisible(trends[, hist(mass.sc, main = 'log Mass (g)', cex.main = cexmain)])
invisible(trends[, hist(speed.sc, main = 'log Speed (km/hr)', cex.main = cexmain)])

invisible(trends[, hist(lifespan.sc, main = 'log Lifespan (yr)', cex.main = cexmain)])
invisible(trends[, hist(consumerfrac.sc, main = 'Consumers (fraction)', cex.main = cexmain)])

invisible(trends[, hist(endothermfrac.sc, main = 'Endotherms (fraction)', cex.main = cexmain)])
invisible(trends[, hist(nspp.sc, main = 'log Species richness', cex.main = cexmain)])

invisible(trends[, hist(thermal_bias.sc, main = 'Thermal bias (°C)', cex.main = cexmain)])
invisible(trends[, hist(npp.sc, main = 'log Net primary productivity', cex.main = cexmain)])

invisible(trends[, hist(human.sc, main = 'Human impact score', cex.main = cexmain)])

Check correlations among variables. Pearson’s r is on the lower diagonal.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- cor(x, y, use = 'pairwise.complete.obs')
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt) #, cex = cex.cor * r)
}
pairs(formula = ~ REALM + tempave.sc + tempave_metab.sc + seas.sc + microclim.sc + temptrend.sc + temptrend_abs.sc +  mass.sc + speed.sc + lifespan.sc + consumerfrac.sc + endothermfrac.sc + nspp.sc + thermal_bias.sc + npp.sc + human.sc, data = trends, gap = 1/10, cex = 0.2, col = '#00000022', lower.panel = panel.cor)

Mass and lifespan look tightly correlated, but r only 0.56…? Tempave_metab and lifespan don’t look tightly correlated, but r= -0.81 Tempave_metab and speed don’t look tightly correlated, but r= -0.83 Lifespan and speed don’t look tightly correlated, but r = 0.73

Examine how many data points are available

Just turnover

cat('Overall # time-series: ', nrow(trends), '\n')

Overall # time-series:  53013 

cat('# studies: ', trends[, length(unique(STUDY_ID))], '\n')

# studies:  332 

cat('Data points: ', trends[, sum(nyrBT)], '\n')

Data points:  293973 

trends[, table(REALM)]

REALM
 Freshwater      Marine Terrestrial 
       1025       48647        3341 

trends[, table(taxa_mod)]

taxa_mod
                        All                  Amphibians                     Benthos                       Birds 
                       1705                         379                        4679                       13741 
                       Fish               Invertebrates                     Mammals Marine invertebrates/plants 
                      28473                        2996                         525                         206 
                      Plant                    Reptiles 
                        305                           4 

trends[, table(taxa_mod, REALM)]

                             REALM
taxa_mod                      Freshwater Marine Terrestrial
  All                                  0   1702           3
  Amphibians                           2      0         377
  Benthos                              0   4679           0
  Birds                                0  11099        2642
  Fish                              1006  27467           0
  Invertebrates                       15   2901          80
  Mammals                              0    478          47
  Marine invertebrates/plants          0    206           0
  Plant                                1    115         189
  Reptiles                             1      0           3

With all covariates

# the cases we can compare
apply(trends[, .(Jtutrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc, temptrend.sc, mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)], MARGIN = 2, FUN = function(x) sum(!is.na(x)))

        Jtutrend            REALM       tempave.sc tempave_metab.sc          seas.sc     microclim.sc 
           53013            53013            49916            49916            49916            51834 
    temptrend.sc          mass.sc         speed.sc      lifespan.sc  consumerfrac.sc endothermfrac.sc 
           49916            52820            52689            51540            47534            53013 
         nspp.sc  thermal_bias.sc           npp.sc         human.sc 
           53013            49371            52863            53013 

i <- trends[, complete.cases(Jtutrend, temptrend.sc, tempave_metab.sc, REALM, seas.sc, microclim.sc, npp.sc, mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, thermal_bias.sc)]
cat('Overall # time-series: ', sum(i), '\n')

Overall # time-series:  43585 

cat('# studies: ', trends[i, length(unique(STUDY_ID))], '\n')

# studies:  250 

cat('Data points: ', trends[i, sum(nyrBT)], '\n')

Data points:  222824 

trends[i, table(REALM)]

REALM
 Freshwater      Marine Terrestrial 
       1008       39735        2842 

trends[i, table(taxa_mod)]

taxa_mod
          All    Amphibians       Benthos         Birds          Fish Invertebrates       Mammals         Plant 
          521            12           590         11803         27372          2567           518           200 
     Reptiles 
            2 

trends[i, table(taxa_mod, REALM)]

               REALM
taxa_mod        Freshwater Marine Terrestrial
  All                    0    520           1
  Amphibians             2      0          10
  Benthos                0    590           0
  Birds                  0   9221        2582
  Fish                 993  26379           0
  Invertebrates         12   2484          71
  Mammals                0    477          41
  Plant                  1     64         135
  Reptiles               0      0           2

Choose the variance structure for mixed effects modles

Try combinations of

• variance scaled to a power of the number of years in the community time-series
• variance scaled to a power of the abs temperature trend
• random intercept for taxa_mod
• random intercept for STUDY_ID
• random slope (abs temperature trend) for taxa_mod
• random slope (abs temperature trend) for STUDY_ID
• random intercept for rarefyID (for overdispersion)

And choose the one with lowest AIC (not run: takes a long time)

# fit models for variance structure
fixed <- formula(Jtutrend ~ REALM + tempave_metab.sc + seas.sc + microclim.sc + npp.sc + temptrend_abs.sc +
                     mass.sc + speed.sc + lifespan.sc + consumerfrac.sc + thermal_bias.sc)
i <- trends[, complete.cases(Jtutrend, REALM, tempave_metab.sc, seas.sc, microclim.sc, npp.sc, temptrend_abs.sc,
                             mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, thermal_bias.sc)]
mods <- vector('list', 0)
mods[[1]] <- gls(fixed, data = trends[i,])
mods[[2]] <- gls(fixed, data = trends[i,], weights = varPower(-0.5, ~nyrBT))
mods[[3]] <- gls(fixed, data = trends[i,], weights = varPower(0.5, ~ abs(temptrend)))

mods[[4]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, control = lmeControl(opt = "optim"))
mods[[5]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, control = lmeControl(opt = "optim"))
mods[[6]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID, control = lmeControl(opt = "optim"))
mods[[7]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, control = lmeControl(opt = "optim"))
mods[[8]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, control = lmeControl(opt = "optim"))

mods[[9]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | taxa_mod)
mods[[10]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | STUDY_ID)
mods[[11]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | taxa_mod2/STUDY_ID, control = lmeControl(opt = "optim"))
mods[[12]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)) # includes overdispersion. new formula so that random slope is only for study level (not enough data to extend to rarefyID).
mods[[13]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)) # 30+ min to fit

mods[[14]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[15]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, weights = varPower(-0.5, ~nyrBT))
mods[[16]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[17]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, weights = varPower(-0.5, ~nyrBT))
mods[[18]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, weights = varPower(-0.5, ~nyrBT))
mods[[19]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[20]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT))
mods[[21]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ 1), weights = varPower(-0.5, ~nyrBT))
mods[[22]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ 1, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT))
mods[[23]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc), weights = varPower(-0.5, ~nyrBT)) # singular precision warning with lmeControl(opt = 'optim') and convergence error without
mods[[24]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT)) # singular precision warning with lmeControl(opt = 'optim') and convergence error without

mods[[25]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, weights = varPower(-0.5, ~abs(temptrend)))
mods[[26]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[27]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[28]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[29]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[30]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|taxa_mod2/STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)), control = lmeControl(opt = "optim"))
mods[[31]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~abs(temptrend)))
mods[[32]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~abs(temptrend)), control = lmeControl(opt = "optim")) # singular precision warning

aics <- sapply(mods, AIC)
minaics <- aics - min(aics)
minaics
which.min(aics)

Chooses the random slopes (temptrend_abs) & intercepts for STUDY_ID, overdispersion, and variance scaled to number of years. We haven’t dealt with potential testing on the boundary issues here yet.

Results

Where do we have data?

world <- map_data('world')
ggplot(world, aes(x = long, y = lat, group = group)) +
    geom_polygon(fill = 'lightgray', color = 'white') +
    geom_point(data = trends, aes(rarefyID_x, rarefyID_y, group = REALM, color = REALM), size = 0.5, alpha = 0.4)  +
    scale_color_brewer(palette="Set1", name = 'Realm') +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        legend.key=element_blank(),
        axis.text=element_text(size=16),
        axis.title=element_text(size=20)) +
  labs(x = 'Longitude (°)', y = 'Latitude (°)')

Mostly northern hemisphere, but spread all over. No so much in Africa or much of Asia.

Plot turnover vs. explanatory variables

Lines are ggplot smoother fits by realm.

Strong trends with temperature change, but trends are pretty symmetric around no trend in temperature, which implies warming or cooling drives similar degree of community turnover. Some indication of less turnover for larger organisms (mass) Higher turnover on land with higher seasonality? More turnover for shorter-lived organisms? No really clear differences among realms.

Nicer plots of turnover vs. temperature?

Violin plots

trends[temptrend <= -0.5, temptrendtext := 'Cooling']
trends[abs(temptrend) <= 0.1, temptrendtext := 'Stable']
trends[temptrend >= 0.5, temptrendtext := 'Warming']

trends[abs(rarefyID_y) < 35, latzone := 'Subtropics']
trends[abs(rarefyID_y) >= 35 & abs(rarefyID_x) < 66.56339, latzone := 'Temperate'] 
trends[abs(rarefyID_y) >= 66.56339, latzone := 'Polar']
trends[, latzone := factor(latzone, levels = c('Subtropics', 'Temperate', 'Polar'))]

p1 <- ggplot(trends, aes(temptrend, Jtutrend, color = REALM, fill = REALM, size = nyrBT)) +
  geom_point(na.rm = TRUE, shape = 16, alpha = 0.1) + 
  geom_smooth(data=subset(trends, abs(temptrend) < 0.75), method = 'gam', formula = y ~ s(x, bs = "cs"), 
              na.rm = TRUE) +
  scale_color_brewer(palette="Set1", name = 'Realm') + 
  scale_fill_brewer(palette="Set1", name = 'Realm') + 
  labs(x = 'Temperature trend (°C/yr)', y = 'Jaccard turnover') +
  scale_size_continuous(range = c(1, 8), breaks = c(2, 5, 20)) +
  guides(size = guide_legend(title = 'Years',
                             override.aes = list(linetype=0, fill = NA, alpha = 1))) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        legend.key=element_blank(),
        legend.key.size = unit(0.5,"line"), 
        axis.text=element_text(size=8),
        axis.title=element_text(size=10))

p2 <- ggplot(trends[!is.na(temptrendtext), ], aes(temptrendtext, Jtutrend)) +
  geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), fill = 'grey') +
  labs(x = '', y = 'Jaccard turnover') +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        legend.key=element_blank(),
        axis.text=element_text(size=8),
        axis.title=element_text(size=10))

p3 <- ggplot(trends[abs(temptrend) >= 0.5 & !is.na(latzone), ], aes(latzone, Jtutrend)) +
  geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), fill = 'grey') + 
  labs(x = '', y = '') +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        legend.key=element_blank(),
        axis.text=element_text(size=7),
        axis.title=element_text(size=10))

grid.arrange(p1, p2, p3, ncol = 2, layout_matrix = rbind(c(1,1), c(2,3)),
             heights=c(unit(0.66, "npc"), unit(0.34, "npc")))

Average rates of turnover

trends[abs(temptrend) >= 0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # turnover per year for locations changing temperature
trends[abs(temptrend) < 0.1, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # not changing temperature
trends[temptrend >= 0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # warming
trends[temptrend <= -0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # cooling

trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) < 35, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # tropics and sub-tropics
trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) >= 35 & abs(rarefyID_y) < 66.56339, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # temperate
trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) >= 66.56339, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # arctic

Hexagonal bins

require(ggplot2)
ggplot(trends[REALM == 'Terrestrial',], aes(temptrend, Jtutrend)) +
  geom_hex() + 
  scale_fill_gradient(trans = "log")

ggplot(trends[REALM == 'Marine',], aes(temptrend, Jtutrend)) +
  geom_hex() + 
  scale_fill_gradient(trans = "log")

ggplot(trends[REALM == 'Freshwater',], aes(temptrend, Jtutrend)) +
  geom_hex() + 
  scale_fill_gradient(trans = "log")

Compare covariates across realms

i <- trends[, !duplicated(rarefyID)]; sum(i)
par(mfrow=c(5,3))
beanplot(rarefyID_y ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Latitude (degN)', ll = 0.05)
beanplot(tempave ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Temperature (degC)', ll = 0.05)
beanplot(tempave_metab ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Metabolic Temperature (degC)', ll = 0.05, bw = 'nrd0') # nrd0 bandwidth to calculation gap
beanplot(seas ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Seasonality (degC)', ll = 0.05)
beanplot(microclim ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Microclimates (degC)', ll = 0.05)
beanplot(temptrend ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Temperature trend (degC/yr)', ll = 0.05)
beanplot(mass_mean_weight ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Mass (g)', ll = 0.05, log = 'y')
beanplot(speed_mean_weight +1 ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Speed (km/hr)', ll = 0.05, log = 'y')
beanplot(lifespan_mean_weight ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Lifespan (yr)', ll = 0.05, log = 'y')
#beanplot(consfrac ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Consumers (fraction)', ll = 0.05, log = '') # too sparse
#beanplot(endofrac ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Endotherms (fraction)', ll = 0.05, log = '') # too sparse
beanplot(Nspp ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Number of species', ll = 0.05, log = 'y')
beanplot(thermal_bias ~ REALM, data = trends[i & !is.na(thermal_bias),], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Thermal bias (degC)', ll = 0.05)
beanplot(npp ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'NPP', ll = 0.05)
beanplot(human ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Human impact score', ll = 0.05)

Marine are in generally warmer locations (seawater doesn’t freeze) Marine have much lower seasonality. Marine and freshwater have some very small masses (plankton), but much of dataset is similar to terrestrial. Marine has a lot of slow, crawling organisms, but land has plants. Land also has birds (fast).

Jaccard turnover temporal trend

Try static covariates plus interactions of abs temperature trend with each covariate:

• realm
• environmental temperature
• average metabolic temperature
• seasonality
• microclimates
• NPP
• speed
• mass
• lifespan
• consumer vs. producer
• thermal bias

Except for thermal bias: interact with temperature trend (not abs)

Temperature-only model (Jtutrend, Jbetatrend, Horntrend)

summary(modonlyTtrend)

Linear mixed-effects model fit by REML
 Data: trends[i, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.04460504 (Intr)
temptrend_abs.sc 0.01904043 0.156 

 Formula: ~1 | rarefyID %in% STUDY_ID
        (Intercept)  Residual
StdDev: 0.001713667 0.2946268

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
   power 
-1.23577 
Fixed effects: Jtutrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.175                                
REALMMarine                     -0.928  0.162                         
REALMTerrestrial                -0.928  0.162  0.861                  
abs(temptrend):REALMMarine       0.173 -0.990 -0.164 -0.161           
abs(temptrend):REALMTerrestrial  0.172 -0.984 -0.160 -0.166  0.974    

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-8.46254351 -0.30730773  0.08443652  0.54660915  6.61661896 

Number of Observations: 49916
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   317                  49916 

summary(modonlyTtrendJbeta)

Linear mixed-effects model fit by REML
 Data: trends[i2, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.05659393 (Intr)
temptrend_abs.sc 0.02520050 0.021 

 Formula: ~1 | rarefyID %in% STUDY_ID
         (Intercept)  Residual
StdDev: 1.693577e-06 0.3590015

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
    power 
-1.413558 
Fixed effects: Jbetatrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.160                                
REALMMarine                     -0.928  0.148                         
REALMTerrestrial                -0.931  0.149  0.864                  
abs(temptrend):REALMMarine       0.158 -0.991 -0.150 -0.147           
abs(temptrend):REALMTerrestrial  0.157 -0.986 -0.146 -0.151  0.976    

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-7.8997139 -0.1696297  0.2037637  0.6704994  6.6449583 

Number of Observations: 49916
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   317                  49916 

summary(modonlyTtrendHorn)

Linear mixed-effects model fit by REML
 Data: trends[i3, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.05792109 (Intr)
temptrend_abs.sc 0.02134495 0.134 

 Formula: ~1 | rarefyID %in% STUDY_ID
        (Intercept) Residual
StdDev: 0.007146559 0.306766

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
    power 
-1.217277 
Fixed effects: Horntrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.154                                
REALMMarine                     -0.914  0.141                         
REALMTerrestrial                -0.925  0.143  0.845                  
abs(temptrend):REALMMarine       0.153 -0.989 -0.142 -0.141           
abs(temptrend):REALMTerrestrial  0.152 -0.984 -0.139 -0.146  0.973    

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-6.16836348 -0.32501087  0.06399723  0.53830958  5.96234033 

Number of Observations: 48800
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   276                  48800 

Try simplifying the temp-only models

if(file.exists('temp/modonlyTtrendsimpreml.rds')){
  modonlyTtrendsimpreml <- readRDS('temp/modonlyTtrendsimpreml.rds')
} else {
  require(MASS) # for stepAIC
  modonlyTtrendml <- update(modonlyTtrend, method = 'ML')
  modonlyTtrendsimp <- stepAIC(modonlyTtrendml, direction = 'backward')
  modonlyTtrendsimpreml <- update(modonlyTtrendsimp, method = 'REML')
  saveRDS(modonlyTtrendsimpreml, file = 'temp/modonlyTtrendsimpreml.rds')
}


if(file.exists('temp/modonlyTtrendJbetasimpreml.rds')){
  modonlyTtrendJbetasimpreml <- readRDS('temp/modonlyTtrendJbetasimpreml.rds')
} else {
  require(MASS) # for stepAIC
  modonlyTtrendJbetaml <- update(modonlyTtrendJbeta, method = 'ML')
  modonlyTtrendJbetasimp <- stepAIC(modonlyTtrendJbetaml, direction = 'backward')
  modonlyTtrendJbetasimpreml <- update(modonlyTtrendJbetasimp, method = 'REML')
  saveRDS(modonlyTtrendJbetasimpreml, file = 'temp/modonlyTtrendJbetasimpreml.rds')
}

Start:  AIC=-165090
Jbetatrend ~ abs(temptrend) * REALM

                       Df     AIC
<none>                    -165090
- abs(temptrend):REALM  2 -164091

if(file.exists('temp/modonlyTtrendHornsimpreml.rds')){
  modonlyTtrendHornsimpreml <- readRDS('temp/modonlyTtrendHornsimpreml.rds')
} else {
  require(MASS) # for stepAIC
  modonlyTtrendHornml <- update(modonlyTtrendHorn, method = 'ML')
  modonlyTtrendHornsimp <- stepAIC(modonlyTtrendHornml, direction = 'backward')
  modonlyTtrendHornsimpreml <- update(modonlyTtrendHornsimp, method = 'REML')
  saveRDS(modonlyTtrendHornsimpreml, file = 'temp/modonlyTtrendHornsimpreml.rds')
}

Start:  AIC=-146055.2
Horntrend ~ abs(temptrend) * REALM

                       Df     AIC
<none>                    -146055
- abs(temptrend):REALM  2 -145334

summary(modonlyTtrendsimpreml)

Linear mixed-effects model fit by REML
 Data: trends[i, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.04460504 (Intr)
temptrend_abs.sc 0.01904043 0.156 

 Formula: ~1 | rarefyID %in% STUDY_ID
        (Intercept)  Residual
StdDev: 0.001713667 0.2946268

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
   power 
-1.23577 
Fixed effects: Jtutrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.175                                
REALMMarine                     -0.928  0.162                         
REALMTerrestrial                -0.928  0.162  0.861                  
abs(temptrend):REALMMarine       0.173 -0.990 -0.164 -0.161           
abs(temptrend):REALMTerrestrial  0.172 -0.984 -0.160 -0.166  0.974    

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-8.46254351 -0.30730773  0.08443652  0.54660915  6.61661896 

Number of Observations: 49916
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   317                  49916 

summary(modonlyTtrendJbetasimpreml)

Linear mixed-effects model fit by REML
 Data: trends[i2, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.05659393 (Intr)
temptrend_abs.sc 0.02520050 0.021 

 Formula: ~1 | rarefyID %in% STUDY_ID
         (Intercept)  Residual
StdDev: 1.693577e-06 0.3590015

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
    power 
-1.413558 
Fixed effects: Jbetatrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.160                                
REALMMarine                     -0.928  0.148                         
REALMTerrestrial                -0.931  0.149  0.864                  
abs(temptrend):REALMMarine       0.158 -0.991 -0.150 -0.147           
abs(temptrend):REALMTerrestrial  0.157 -0.986 -0.146 -0.151  0.976    

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-7.8997139 -0.1696297  0.2037637  0.6704994  6.6449583 

Number of Observations: 49916
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   317                  49916 

summary(modonlyTtrendHornsimpreml)

Linear mixed-effects model fit by REML
 Data: trends[i3, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.05792109 (Intr)
temptrend_abs.sc 0.02134495 0.134 

 Formula: ~1 | rarefyID %in% STUDY_ID
        (Intercept) Residual
StdDev: 0.007146559 0.306766

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
    power 
-1.217277 
Fixed effects: Horntrend ~ abs(temptrend) * REALM 

 Correlation: 
                                (Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend)                  -0.154                                
REALMMarine                     -0.914  0.141                         
REALMTerrestrial                -0.925  0.143  0.845                  
abs(temptrend):REALMMarine       0.153 -0.989 -0.142 -0.141           
abs(temptrend):REALMTerrestrial  0.152 -0.984 -0.139 -0.146  0.973    

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-6.16836348 -0.32501087  0.06399723  0.53830958  5.96234033 

Number of Observations: 48800
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   276                  48800 

Plot the temp-only coefficients

coefs <- summary(modonlyTtrend)$tTable
coefs2 <- summary(modonlyTtrendJbeta)$tTable
coefs3 <- summary(modonlyTtrendHorn)$tTable

par(las = 1, mai = c(0.8, 2, 0.1, 0.1))
rows1 <- which(grepl('temptrend', rownames(coefs)))
plot(0,0, col = 'white', xlim=c(-0.02, 0.85), ylim = c(1,length(rows1)+0.3), 
     yaxt='n', xlab = 'Turnover per |°C/yr|', ylab ='')
axis(2, at = length(rows1):1, labels = c('Freshwater', 'Marine', 'Terrestrial'), cex.axis = 0.7)

abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
  x = coefs[rows1[i], 1]
  x2 = coefs2[rows1[i], 1]
  x3 = coefs3[rows1[i], 1]
  if(i>1){ # if an interaction, add the intercept
    x = x + coefs['abs(temptrend)', 1]
    x2 = x2 + coefs2['abs(temptrend)', 1]
    x3 = x3 + coefs3['abs(temptrend)', 1]
  }
  se = coefs[rows1[i], 2]
  se2 = coefs2[rows1[i], 2]
  se3 = coefs3[rows1[i], 2]
  points(x, length(rows1) + 1 - i, pch = 16)
  lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i, length(rows1) + 1 - i))
  
  points(x2, length(rows1) + 1.1 - i, pch = 16, col = 'light grey')
  lines(x = c(x2-se2, x2+se2), y = c(length(rows1) + 1.1 - i, length(rows1) + 1.1 - i), col = 'light grey')
  
  points(x3, length(rows1) + 1.2 - i, pch = 16, col = 'dark grey')
  lines(x = c(x3-se3, x3+se3), y = c(length(rows1) + 1.2 - i, length(rows1) + 1.2 - i), col = 'dark grey')
}

legend('bottomright', col = c('black', 'dark grey', 'light grey'), lwd = 1, pch = 16, 
       legend = c('Jaccard turnover', 'Jaccard total', 'Horn-Morisita'))

Full model

Try simplifying the model

summary(modTfull1simpreml)

Linear mixed-effects model fit by REML
 Data: trends[i, ] 

Random effects:
 Formula: ~temptrend_abs.sc | STUDY_ID
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
(Intercept)      0.05008648 (Intr)
temptrend_abs.sc 0.01884335 0.56  

 Formula: ~1 | rarefyID %in% STUDY_ID
        (Intercept)  Residual
StdDev: 0.001194902 0.3037795

Variance function:
 Structure: Power of variance covariate
 Formula: ~nyrBT 
 Parameter estimates:
    power 
-1.246946 
Fixed effects: Jtutrend ~ temptrend_abs.sc + REALM + tempave.sc + tempave_metab.sc +      seas.sc + microclim.sc + mass.sc + speed.sc + lifespan.sc +      endothermfrac.sc + nspp.sc + temptrend.sc + thermal_bias.sc +      npp.sc + human.sc + temptrend_abs.sc:REALM + temptrend_abs.sc:tempave_metab.sc +      temptrend_abs.sc:seas.sc + temptrend_abs.sc:microclim.sc +      temptrend_abs.sc:mass.sc + temptrend_abs.sc:lifespan.sc +      temptrend_abs.sc:nspp.sc + temptrend.sc:thermal_bias.sc +      temptrend_abs.sc:npp.sc + temptrend_abs.sc:human.sc + REALM:human.sc +      temptrend_abs.sc:REALM:human.sc 

 Correlation: 
                                           (Intr) tmpt_. REALMMr REALMTr tmpv.s tmpv_. ses.sc mcrcl. mss.sc spd.sc
temptrend_abs.sc                            0.360                                                                 
REALMMarine                                -0.918 -0.333                                                          
REALMTerrestrial                           -0.916 -0.325  0.835                                                   
tempave.sc                                  0.004  0.005 -0.014  -0.005                                           
tempave_metab.sc                            0.052  0.023 -0.037  -0.023  -0.289                                   
seas.sc                                    -0.036  0.000  0.053  -0.006   0.268 -0.084                            
microclim.sc                               -0.005  0.007  0.012  -0.016   0.049  0.063  0.171                     
mass.sc                                     0.022  0.000 -0.009  -0.020   0.006 -0.551 -0.045 -0.015              
speed.sc                                    0.013  0.005 -0.025   0.022  -0.003  0.154 -0.072 -0.034 -0.128       
lifespan.sc                                 0.032  0.024 -0.015  -0.006  -0.047  0.750  0.095  0.046 -0.807  0.246
endothermfrac.sc                            0.145  0.007 -0.055  -0.222   0.123 -0.212  0.070  0.017  0.037  0.023
nspp.sc                                    -0.026 -0.009 -0.002  -0.014  -0.042 -0.013 -0.032 -0.053 -0.161 -0.043
temptrend.sc                               -0.006 -0.013  0.006   0.007  -0.050 -0.018 -0.041 -0.014  0.010 -0.009
thermal_bias.sc                             0.035  0.007 -0.047  -0.016   0.760 -0.053 -0.156 -0.135  0.043  0.016
npp.sc                                      0.001 -0.002 -0.003   0.015  -0.322  0.237 -0.255 -0.221 -0.009 -0.032
human.sc                                   -0.167 -0.020  0.157   0.152   0.018 -0.005  0.009  0.017  0.005  0.004
temptrend_abs.sc:REALMMarine               -0.341 -0.955  0.364   0.308   0.000 -0.009  0.011 -0.002  0.004 -0.001
temptrend_abs.sc:REALMTerrestrial          -0.319 -0.896  0.299   0.337  -0.029 -0.021 -0.027 -0.024  0.008 -0.001
temptrend_abs.sc:tempave_metab.sc           0.025  0.081 -0.008  -0.016   0.103  0.365  0.049  0.080 -0.302  0.047
temptrend_abs.sc:seas.sc                    0.019 -0.014 -0.010  -0.039   0.126  0.026  0.222  0.029 -0.029 -0.017
temptrend_abs.sc:microclim.sc              -0.005 -0.036  0.001  -0.001   0.017  0.100  0.060  0.312 -0.036  0.015
temptrend_abs.sc:mass.sc                    0.002 -0.001  0.003  -0.004  -0.008 -0.273 -0.037 -0.025  0.457 -0.030
temptrend_abs.sc:lifespan.sc                0.009  0.050  0.001   0.003  -0.011  0.367  0.050  0.022 -0.385  0.064
temptrend_abs.sc:nspp.sc                   -0.008  0.021 -0.002   0.001  -0.028  0.014  0.001  0.013 -0.070  0.029
temptrend.sc:thermal_bias.sc                0.009  0.019 -0.008  -0.007   0.067  0.033  0.000 -0.032 -0.002  0.008
temptrend_abs.sc:npp.sc                     0.002 -0.054 -0.007   0.012  -0.115  0.095 -0.161 -0.163  0.013  0.000
temptrend_abs.sc:human.sc                  -0.070  0.160  0.066   0.063  -0.004 -0.009  0.005  0.006  0.009 -0.007
REALMMarine:human.sc                        0.167  0.021 -0.158  -0.152  -0.020  0.004 -0.017 -0.024 -0.005 -0.003
REALMTerrestrial:human.sc                   0.166  0.020 -0.157  -0.152  -0.032  0.007 -0.022 -0.006 -0.002 -0.001
temptrend_abs.sc:REALMMarine:human.sc       0.070 -0.158 -0.066  -0.063   0.001  0.009 -0.009 -0.009 -0.009  0.007
temptrend_abs.sc:REALMTerrestrial:human.sc  0.069 -0.160 -0.066  -0.063   0.012  0.007 -0.003 -0.005 -0.008  0.006
                                           lfspn. endth. nspp.s tmptr. thrm_. npp.sc hmn.sc tm_.:REALMM tm_.:REALMT
temptrend_abs.sc                                                                                                   
REALMMarine                                                                                                        
REALMTerrestrial                                                                                                   
tempave.sc                                                                                                         
tempave_metab.sc                                                                                                   
seas.sc                                                                                                            
microclim.sc                                                                                                       
mass.sc                                                                                                            
speed.sc                                                                                                           
lifespan.sc                                                                                                        
endothermfrac.sc                           -0.054                                                                  
nspp.sc                                     0.107  0.051                                                           
temptrend.sc                               -0.007 -0.002  0.013                                                    
thermal_bias.sc                            -0.075  0.028 -0.100  0.003                                             
npp.sc                                      0.053 -0.129 -0.161 -0.013 -0.139                                      
human.sc                                    0.008  0.006 -0.015  0.008  0.015 -0.029                               
temptrend_abs.sc:REALMMarine               -0.011  0.005 -0.003  0.008 -0.003 -0.006  0.020                        
temptrend_abs.sc:REALMTerrestrial          -0.017  0.019  0.000  0.015 -0.023  0.020  0.019  0.855                 
temptrend_abs.sc:tempave_metab.sc           0.385  0.027 -0.011 -0.049  0.108  0.003 -0.001 -0.048      -0.098     
temptrend_abs.sc:seas.sc                    0.035  0.053  0.001 -0.011  0.087 -0.178 -0.003  0.046      -0.082     
temptrend_abs.sc:microclim.sc               0.025 -0.043  0.003  0.000 -0.026 -0.142  0.002  0.044       0.002     
temptrend_abs.sc:mass.sc                   -0.381  0.023 -0.060  0.013  0.017  0.008  0.010  0.011       0.017     
temptrend_abs.sc:lifespan.sc                0.466 -0.031  0.053 -0.010 -0.027  0.011 -0.004 -0.032      -0.038     
temptrend_abs.sc:nspp.sc                    0.076 -0.029  0.303  0.015 -0.051 -0.017 -0.011 -0.053      -0.039     
temptrend.sc:thermal_bias.sc                0.003  0.008 -0.010 -0.402 -0.015 -0.003 -0.006 -0.017      -0.020     
temptrend_abs.sc:npp.sc                     0.010 -0.062 -0.034 -0.016 -0.044  0.366 -0.012  0.040       0.053     
temptrend_abs.sc:human.sc                  -0.008  0.000 -0.003  0.007 -0.008 -0.012  0.409 -0.151      -0.139     
REALMMarine:human.sc                       -0.009 -0.006  0.013 -0.008 -0.014  0.018 -0.998 -0.021      -0.019     
REALMTerrestrial:human.sc                  -0.013 -0.011  0.018 -0.007 -0.020  0.027 -0.998 -0.020      -0.019     
temptrend_abs.sc:REALMMarine:human.sc       0.008  0.000  0.000 -0.006  0.006  0.008 -0.408  0.150       0.139     
temptrend_abs.sc:REALMTerrestrial:human.sc  0.008  0.002  0.004 -0.007  0.011  0.013 -0.406  0.151       0.138     
                                           t_.:_. tmptrnd_bs.sc:s. tmptrnd_bs.sc:mc. tmptrnd_bs.sc:ms.
temptrend_abs.sc                                                                                      
REALMMarine                                                                                           
REALMTerrestrial                                                                                      
tempave.sc                                                                                            
tempave_metab.sc                                                                                      
seas.sc                                                                                               
microclim.sc                                                                                          
mass.sc                                                                                               
speed.sc                                                                                              
lifespan.sc                                                                                           
endothermfrac.sc                                                                                      
nspp.sc                                                                                               
temptrend.sc                                                                                          
thermal_bias.sc                                                                                       
npp.sc                                                                                                
human.sc                                                                                              
temptrend_abs.sc:REALMMarine                                                                          
temptrend_abs.sc:REALMTerrestrial                                                                     
temptrend_abs.sc:tempave_metab.sc                                                                     
temptrend_abs.sc:seas.sc                    0.170                                                     
temptrend_abs.sc:microclim.sc              -0.012  0.058                                              
temptrend_abs.sc:mass.sc                   -0.626 -0.016            0.009                             
temptrend_abs.sc:lifespan.sc                0.829  0.045            0.021            -0.813           
temptrend_abs.sc:nspp.sc                    0.080 -0.056           -0.174            -0.224           
temptrend.sc:thermal_bias.sc                0.072 -0.030           -0.036            -0.004           
temptrend_abs.sc:npp.sc                     0.046 -0.169           -0.213            -0.032           
temptrend_abs.sc:human.sc                  -0.030 -0.018            0.038             0.049           
REALMMarine:human.sc                        0.001 -0.001           -0.005            -0.009           
REALMTerrestrial:human.sc                  -0.001 -0.002            0.001            -0.008           
temptrend_abs.sc:REALMMarine:human.sc       0.030  0.009           -0.037            -0.048           
temptrend_abs.sc:REALMTerrestrial:human.sc  0.023  0.014           -0.027            -0.045           
                                           tmptrnd_bs.sc:l. tmptrnd_bs.sc:ns. tm.:_. tmptrnd_bs.sc:np.
temptrend_abs.sc                                                                                      
REALMMarine                                                                                           
REALMTerrestrial                                                                                      
tempave.sc                                                                                            
tempave_metab.sc                                                                                      
seas.sc                                                                                               
microclim.sc                                                                                          
mass.sc                                                                                               
speed.sc                                                                                              
lifespan.sc                                                                                           
endothermfrac.sc                                                                                      
nspp.sc                                                                                               
temptrend.sc                                                                                          
thermal_bias.sc                                                                                       
npp.sc                                                                                                
human.sc                                                                                              
temptrend_abs.sc:REALMMarine                                                                          
temptrend_abs.sc:REALMTerrestrial                                                                     
temptrend_abs.sc:tempave_metab.sc                                                                     
temptrend_abs.sc:seas.sc                                                                              
temptrend_abs.sc:microclim.sc                                                                         
temptrend_abs.sc:mass.sc                                                                              
temptrend_abs.sc:lifespan.sc                                                                          
temptrend_abs.sc:nspp.sc                    0.231                                                     
temptrend.sc:thermal_bias.sc               -0.004            0.002                                    
temptrend_abs.sc:npp.sc                     0.042           -0.188             0.001                  
temptrend_abs.sc:human.sc                  -0.027           -0.032            -0.004 -0.034           
REALMMarine:human.sc                        0.003            0.010             0.005  0.008           
REALMTerrestrial:human.sc                   0.002            0.013             0.004  0.012           
temptrend_abs.sc:REALMMarine:human.sc       0.027            0.028             0.005  0.024           
temptrend_abs.sc:REALMTerrestrial:human.sc  0.021            0.033             0.006  0.030           
                                           tmptrnd_bs.sc:h. REALMM: REALMT: t_.:REALMM:
temptrend_abs.sc                                                                       
REALMMarine                                                                            
REALMTerrestrial                                                                       
tempave.sc                                                                             
tempave_metab.sc                                                                       
seas.sc                                                                                
microclim.sc                                                                           
mass.sc                                                                                
speed.sc                                                                               
lifespan.sc                                                                            
endothermfrac.sc                                                                       
nspp.sc                                                                                
temptrend.sc                                                                           
thermal_bias.sc                                                                        
npp.sc                                                                                 
human.sc                                                                               
temptrend_abs.sc:REALMMarine                                                           
temptrend_abs.sc:REALMTerrestrial                                                      
temptrend_abs.sc:tempave_metab.sc                                                      
temptrend_abs.sc:seas.sc                                                               
temptrend_abs.sc:microclim.sc                                                          
temptrend_abs.sc:mass.sc                                                               
temptrend_abs.sc:lifespan.sc                                                           
temptrend_abs.sc:nspp.sc                                                               
temptrend.sc:thermal_bias.sc                                                           
temptrend_abs.sc:npp.sc                                                                
temptrend_abs.sc:human.sc                                                              
REALMMarine:human.sc                       -0.408                                      
REALMTerrestrial:human.sc                  -0.408            0.997                     
temptrend_abs.sc:REALMMarine:human.sc      -0.998            0.409   0.407             
temptrend_abs.sc:REALMTerrestrial:human.sc -0.993            0.405   0.404   0.991     

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-7.59387150 -0.34529883  0.05138559  0.55509022  6.91268058 

Number of Observations: 43585
Number of Groups: 
              STUDY_ID rarefyID %in% STUDY_ID 
                   250                  43585 

Plot the coefficients from the simplified model

coefs <- summary(modTfull1simpreml)$tTable
par(las = 1, mai = c(0.5, 3, 0.1, 0.1))
rows1 <- which(!grepl('Intercept', rownames(coefs)))
plot(0,0, col = 'white', xlim=c(-0.02, 0.08), ylim = c(1,length(rows1)), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows1):1, labels = rownames(coefs)[rows1], cex.axis = 0.7)

abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
  x = coefs[rows1[i], 1]
  se = coefs[rows1[i], 2]
  points(x, length(rows1) + 1 - i, pch = 16)
  lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i, length(rows1) + 1 - i))
}

Plot residuals against each predictor

resids <- resid(modTfull1simpreml)
preds <- getData(modTfull1simpreml)
col = '#00000033'
cex = 0.5
par(mfrow = c(4,4))
boxplot(resids ~ preds$REALM, cex = cex, col = col)
plot(preds$temptrend_abs.sc, resids, cex = cex, col = col)
plot(preds$temptrend.sc, resids, cex = cex, col = col)
plot(preds$tempave.sc, resids, cex = cex, col = col)
plot(preds$tempave_metab.sc, resids, cex = cex, col = col)
plot(preds$seas.sc, resids, cex = cex, col = col)
plot(preds$microclim.sc, resids, cex = cex, col = col)
plot(preds$mass.sc, resids, cex = cex, col = col)
plot(preds$speed.sc, resids, cex = cex, col = col)
plot(preds$lifespan.sc, resids, cex = cex, col = col)
plot(preds$consumerfrac.sc, resids, cex = cex, col = col)
plot(preds$endothermfrac.sc, resids, cex = cex, col = col)
plot(preds$nspp.sc, resids, cex = cex, col = col)
plot(preds$thermal_bias.sc, resids, cex = cex, col = col)
plot(preds$npp.sc, resids, cex = cex, col = col)
plot(preds$human.sc, resids, cex = cex, col = col)

Sensitivity analysis: total turnover and Morisita-Horn models

i2 <- trends[, complete.cases(Jbetatrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc, 
                             temptrend.sc, temptrend_abs.sc, mass.sc, speed.sc, lifespan.sc, 
                             consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)]
i3 <- trends[, complete.cases(Horntrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc, 
                             temptrend.sc, temptrend_abs.sc, mass.sc, speed.sc, lifespan.sc, 
                             consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)]

randef <- list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)
varef <- varPower(-0.5, ~nyrBT)

# full models
if(file.exists('temp/modTfullJbeta.rds')){
  modTfullJbeta <- readRDS('temp/modTfullJbeta.rds')
} else {
  modTfullJbeta <- lme(Jbetatrend ~ temptrend_abs.sc*REALM + 
                     temptrend_abs.sc*tempave.sc +
                     temptrend_abs.sc*tempave_metab.sc + 
                     temptrend_abs.sc*seas.sc + 
                     temptrend_abs.sc*microclim.sc + 
                     temptrend_abs.sc*mass.sc + 
                     temptrend_abs.sc*speed.sc + 
                     temptrend_abs.sc*lifespan.sc + 
                     temptrend_abs.sc*consumerfrac.sc +
                     temptrend_abs.sc*endothermfrac.sc +
                     temptrend_abs.sc*nspp.sc +
                     temptrend.sc*thermal_bias.sc +
                     temptrend_abs.sc*npp.sc +
                     temptrend_abs.sc*human.sc*REALM,
                   random = randef, weights = varef, data = trends[i2,], method = 'REML')
  saveRDS(modTfullJbeta, file = 'temp/modTfullJbeta.rds')
}

if(file.exists('temp/modTfullHorn.rds')){
  modTfullHorn <- readRDS('temp/modTfullHorn.rds')
} else {
  modTfullHorn <- lme(Horntrend ~ temptrend_abs.sc*REALM + 
                     temptrend_abs.sc*tempave.sc +
                     temptrend_abs.sc*tempave_metab.sc + 
                     temptrend_abs.sc*seas.sc + 
                     temptrend_abs.sc*microclim.sc + 
                     temptrend_abs.sc*mass.sc + 
                     temptrend_abs.sc*speed.sc + 
                     temptrend_abs.sc*lifespan.sc + 
                     temptrend_abs.sc*consumerfrac.sc +
                     temptrend_abs.sc*endothermfrac.sc +
                     temptrend_abs.sc*nspp.sc +
                     temptrend.sc*thermal_bias.sc +
                     temptrend_abs.sc*npp.sc +
                     temptrend_abs.sc*human.sc*REALM,
                   random = randef, weights = varef, data = trends[i3,], method = 'REML')
  saveRDS(modTfullHorn, file = 'temp/modTfullHorn.rds')
}

summary(modTfullJbeta)
summary(modTfullHorn)

# simplify
if(file.exists('temp/modTfullJbetasimp.rds')){
  modTfullJbetasimp <- readRDS('temp/modTfullJbetasimp.rds')
} else {
  require(MASS) # for stepAIC
  modTfullJbetaml <- update(modTfullJbeta, method = 'ML')
  modTfullJbetasimp <- stepAIC(modTfullJbetaml, direction = 'backward')
  saveRDS(modTfullJbetasimp, file = 'temp/modTfullJbetasimp.rds')
}

if(file.exists('temp/modTfullHornsimp.rds')){
  modTfullHornsimp <- readRDS('temp/modTfullHornsimp.rds')
} else {
  require(MASS) # for stepAIC
  modTfullHornml <- update(modTfullHorn, method = 'ML')
  modTfullHornsimp <- stepAIC(modTfullHornml, direction = 'backward')
  saveRDS(modTfullHornsimp, file = 'temp/modTfullHornsimp.rds')
}

summary(modTfullJbetasimp)
summary(modTfullHornsimp)

if(file.exists('temp/modTfullJbetasimpreml.rds')){
  modTfullJbetasimpreml <- readRDS('temp/modTfullJbetasimpreml.rds')
} else {
  modTfullJbetasimpreml <- update(modTfullJbetasimp, method = 'REML')
  saveRDS(modTfullJbetasimpreml, file = 'temp/modTfullJbetasimpreml.rds')
}
if(file.exists('temp/modTfullHornsimpreml.rds')){
  modTfullHornsimpreml <- readRDS('temp/modTfullHornsimpreml.rds')
} else {
  modTfullHornsimpreml <- update(modTfullHornsimp, method = 'REML')
  saveRDS(modTfullHornsimpreml, file = 'temp/modTfullHornsimpreml.rds')
}

# plot coefs
coefs2 <- summary(modTfullJbetasimpreml)$tTable
coefs3 <- summary(modTfullHornsimpreml)$tTable
varstoplot <- unique(c(rownames(coefs2), rownames(coefs3)))

rows1 <- which(!grepl('Intercept|REALM', varstoplot) | grepl(':', varstoplot)) # vars to plot in first graph
rows1_2 <- which(rownames(coefs2) %in% varstoplot[rows1]) # rows in coefs2
rows1_3 <- which(rownames(coefs3) %in% varstoplot[rows1]) # rows in coefs3
xlims1 <- range(c(coefs2[rows1_2,1] - coefs2[rows1_2,2], 
                  coefs2[rows1_2,1] + coefs2[rows1_2,2], 
                  coefs3[rows1_3,1] - coefs3[rows1_3,2], 
                  coefs3[rows1_3,1] + coefs3[rows1_3,2]))

rows2 <- which(grepl('REALM', varstoplot) & !grepl(':', varstoplot)) # vars to plot in 2nd graph
rows2_2 <- which(rownames(coefs2) %in% varstoplot[rows2]) # rows in coefs2
rows2_3 <- which(rownames(coefs3) %in% varstoplot[rows2]) # rows in coefs3
xlims2 <- range(c(coefs2[rows2_2,1] - coefs2[rows2_2,2], 
                  coefs2[rows2_2,1] + coefs2[rows2_2,2], 
                  coefs3[rows2_3,1] - coefs3[rows2_3,2], 
                  coefs3[rows2_3,1] + coefs3[rows2_3,2]))

cols <- c('black', 'grey') # for Jbeta and Horn models, respectively
offs1 <- 0.1 # offset vertically for the two models
offs2 <- 0.01 # offset vertically for the two models (plot 2)

par(mfrow=c(1,2), las = 1, mai = c(0.5, 3, 0.1, 0.1))

plot(0,0, col = 'white', xlim=xlims1, ylim = c(1,length(rows1)), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows1):1, labels = varstoplot[rows1], cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
  if(varstoplot[rows1[i]] %in% rownames(coefs2)){
    x = coefs2[rownames(coefs2) == varstoplot[rows1[i]], 1]
    se = coefs2[rownames(coefs2) == varstoplot[rows1[i]], 2]
    points(x, length(rows1) + 1 - i + offs1, pch = 16, col = cols[1])
    lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i + offs1, length(rows1) + 1 - i + offs1), col = cols[1])
  }
  if(varstoplot[rows1[i]] %in% rownames(coefs3)){
    x = coefs3[rownames(coefs3) == varstoplot[rows1[i]], 1]
    se = coefs3[rownames(coefs3) == varstoplot[rows1[i]], 2]
    points(x, length(rows1) + 1 - i - offs1, pch = 16, col = cols[2])
    lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i - offs1, length(rows1) + 1 - i - offs1), col = cols[2])
  }
}

plot(0,0, col = 'white', xlim=xlims2, ylim = c(0.9, length(rows2) + 0.1), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows2):1, labels = varstoplot[rows2], cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows2)){
  if(varstoplot[rows2[i]] %in% rownames(coefs2)){
    x = coefs2[rownames(coefs2) == varstoplot[rows2[i]], 1]
    se = coefs2[rownames(coefs2) == varstoplot[rows2[i]], 2]
    points(x, length(rows2) + 1 - i + offs2, pch = 16, col = cols[1])
    lines(x = c(x-se, x+se), y = c(length(rows2) + 1 - i + offs2, length(rows2) + 1 - i + offs2), col = cols[1])
  }
  if(varstoplot[rows2[i]] %in% rownames(coefs3)){
    x = coefs3[rownames(coefs3) == varstoplot[rows2[i]], 1]
    se = coefs3[rownames(coefs3) == varstoplot[rows2[i]], 2]
    points(x, length(rows2) + 1 - i - offs2, pch = 16, col = cols[2])
    lines(x = c(x-se, x+se), y = c(length(rows2) + 1 - i - offs2, length(rows2) + 1 - i - offs2), col = cols[2])
  }
}

Black is for Jaccard total turnover (pres/abs), grey is for Morisita-Horn turnover (considers abundance)

To do